function varargout = normprc (P, pin)

   ## usage:  varargout = normprc (P, pin)
   ##
   ## normalize precipitation series P

   global PAR

   Inanin = isnan(P) ;
   Iinfin = isinf(P) ;

### zscore via NQT ;; FIXME: check Bogner et al. (2012)
   [nr nc] = size(P) ;

   if (nargin < 2)

      CJ = mat2cell(P, nr, ones(1, nc)) ;
      J = num2cell(1:nc) ;
      wrn = warning("query", "Octave:divide-by-zero") ; warning("off", "Octave:divide-by-zero") ;
      [pout Z] = parcellfun_(@(x,j) _normprc({x}, j), CJ, J, PAR.par_opt{:}) ;
      warning(wrn.state, "Octave:divide-by-zero") ;
      Z = cell2mat(Z) ;

      Inanout = isnan(Z) ;
      Iinfout = isinf(Z) ;

      if any((I = (Inanin != Inanout))(:))
	 [i, j] = find(I) ;
	 error("normprc: NaNs created at i, j = %d, %d.\n", i(1), j(1)) ;
      endif 
      if any((I = (Iinfin != Iinfout))(:))
      	 [i, j] = find(I) ;
      	 for l=1:length(i)
      	    warning("xds:xds", "normprc: Infs created at i, j = %d, %d.\n", i(l), j(l)) ;
      	 endfor 
      	 Z(I) = norminv(1 - 1/rows(Z)) ; # replace inf with maximum normal
      endif

      varargout{1} = cell2mat(pout) ;
      varargout{2} = Z ;

   else

      CJ = {} ; for i=1:nc, CJ = {CJ{:}, {P(:,i), pin(:,i)}} ; endfor 
      pout = parcellfun_(@(x) _normprc(x), CJ, PAR.par_opt{:}) ;
      Z = cell2mat(pout) ;

      Inanout = isnan(Z) ;
      Iinfout = isinf(Z) ;

      if (any(any(I = (Inanin != Inanout))))
	 [i, j] = find(I) ;
	 error("normprc: NaNs created at i, j = %d, %d.\n", i(1), j(1)) ;
      endif 
      if (any(any(I = (Iinfin != Iinfout))))
	 [i, j] = find(I) ;
	 for l=1:length(i)
	    warning("xds:xds", "normprc: Infs created at i, j = %d, %d.\n", i(l), j(l)) ;
	 endfor 
	 Z(I) = norminv(1 - 1/rows(Z)) ; # replace inf with maximum normal
      endif 

      varargout{1} = Z ;

   endif 

endfunction


function [res z] = _normprc (C, j)

   ## usage:  [res z] = _normprc (C, j)
   ##
   ## wrapper for normprc

   persistent nfun = 0 ;
   global PAR

   pkg load nan

   nfun++ ;

   P = C{1} ;
   if length(C) > 1, pin = C{2} ; endif 

   I = !isnan(P) ;
   if !any(I)
      error("_normprc: all values undefined.\n") ;
      res = nan ; z = nan(size(P)) ;
      return ;
   endif
   Pm = repmat(mean(P)*PAR.Pzero, rows(P), 1) ; II = (P > Pm) ;
   if any(II)
      x = P(II) ;
      q = sum(I & !II) ./ sum(I) ;
      x0 = min(x) ;
   else
      x = x0 = NaN ; q = 0 ;
   endif

   [F, IF] = cdf(x) ;
   F = (1-q) * F + q ;
   F = norminv(F) ;
   X = [x(IF), F(IF)] ;

   if length(C) == 1

      res = [q ; x0] ;

      if strcmp(PAR.cdffit, "wbl")
	 p0 = [mean(X(:,1)) 0.5 0.5] ;
      else
	 p0 = [0.01 mean(X(:,1)) std(X(:,1))] ;
      endif

      if any(II)

	 opt = optimset("TolX", 1e-10, "TolFun", 1e-10, "MaxFunEvals", 1000, "AutoScaling", "on") ;

	 ##[p, fy, cvg, outp] = nonlin_curvefit(@wblfit, p0, X(:,1), X(:,2)) ;
	 [p v flag out] = fminunc (@(p) cost(p, PAR.cdffit, X, q, x0, PAR.prbwgt), p0, opt) ;
	 [v s1] = cost(p, PAR.cdffit, X, q, x0, PAR.prbwgt) ;

	 if (flag != 1 || s1 > 0.3)
	    [p v flag out] = fminunc (@(p) cost(p, PAR.cdffit, X, q, x0, PAR.prbwgt), p, opt) ;
	    [v s1] = cost(p, PAR.cdffit, X, q, x0, PAR.prbwgt) ;
	    if (flag != 1 && s1 > 0.3) # try harder
	       pkg load optim
	       warning("xds:xds", "normprc: fminunc: (flag, fmin, rms) = (%d, %f, %f)", flag, v, s1) ;
	       opls.ftol = optimget(opt, "TolFun") ;
	       opls.maxev = optimget(opt, "MaxFunEvals") ;
	       [pw v nev] = minimize(@cost, {p, PAR.cdffit, X, q, x0, PAR.prbwgt}, opls) ;
	       [v s1w] = cost(pw, PAR.cdffit, X, q, x0, PAR.prbwgt) ;
	       if s1w > 0.3
		  [pw v nev] = minimize(@cost, {pw, PAR.cdffit, X, q, x0, PAR.prbwgt}, opls) ;
		  [v s1w] = cost(p, PAR.cdffit, X, q, x0, PAR.prbwgt) ;
		  if s1w > 0.3
		     save -binary err.ob p X q x0 PAR opls PAR
		     prb_print(X, p, q, x0) ;
		     warning("normprc (%d): minimize(%d): cost = %f,\nsaving to err.ob\n", j, nfun, s1) ;
		  endif
	       endif
	       pkg unload optim
	       if s1w < s1, p = pw ; endif
	    endif
	 endif

	 if PAR.dbg, prb_print(X, p, q, x0) ; endif

      else

	 p = nan(size(p)) ;

      endif 

      res = [res; p(:)] ;

   else

      q = pin(1) ; x0 = pin(2) ; p = pin(3:end) ;

   endif 

   z = nan(rows(P),1) ;
   if any(II)
      if PAR.prbraw
	 z(II) = F ;
      else
	 z(II) = prb(PAR.cdffit, x, p, q, x0) ;
      endif 
      if (!isreal(z(II)))
	 error("normprc: complex zscores.\n") ;
      endif 
   endif 

   z(I & !II) = trunc_rnd("norm", sum((I & !II)(:)), -Inf, norminv(q)) ;

   if (length(C) > 1)
      res = z ;
   endif 

endfunction


function prb_print (X, p, q, x0, pfile)

   ## usage:  prb_print (X, p, q, x0, pfile)
   ##
   ## print probit estimate

   persistent nprint = 0 ;
   global PAR PRJ BATCH ;
   nprint++ ;

   if (nargin < 5) 
      pdir = fullfile("data", PRJ, PAR.pdd) ;
      pfile = fullfile(pdir, ["prc_" num2str(nprint) ".png"]) ;
   endif

   if BATCH
      graphics_toolkit gnuplot ; figure("visible", "off") ;
   endif
   subplot(2,1,1) ;
   line(X(:,1), X(:,2), "linestyle", "none", "marker", "o", "markersize", 4) ;
   line(X(:,1), feval("prb", PAR.cdffit, X(:,1), p, q, x0), "linestyle", "-", "marker", "none") ;
   xlabel("phys.") ; ylabel("norm.") ;
   subplot(2,1,2) ;
   line(X(:,2), X(:,1), "linestyle", "none", "marker", "o", "markersize", 4) ;
   line(X(:,2), feval("prb_inv", PAR.cdffit, X(:,2), p, q, x0), "linestyle", "-", "marker", "none") ;
   xlabel("norm.") ; ylabel("phys.") ;
   if exist(pdir, "dir"), print(pfile) ; endif
   close ;

   printf("prb_print: --> %s\n", pfile) ;

endfunction
